Molecular Synchronization Waves in Arrays of Allosterically Regulated Enzymes 
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Spatiotemporal pattern formation in a product-activated enzymic reaction at high enzyme con- 
centrations is investigated. Stochastic simulations show that catalytic turnover cycles of individual 
enzymes can become coherent and that complex wave patterns of molecular synchronization can 
develop. The analysis based on the mean-field approximation indicates that the observed patterns 
result from the presence of Hopf and wave bifurcations in the considered system. 
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Molecular machines, such as molecular motors, ion 
pumps and some enzymes, play a fundamental role in 
biological cells and can be also used in the emerging soft- 
matter nanotechnology |l| . A protein machine is a cyclic 
device, where each cycle consists of conformational mo- 
tions initiated by binding of an energy-bringing ligand 
0, 0). In motors, such internal motions generate me- 
chanical work while in enzymes they enable or facil- 
itate chemical reaction events (see, e.g., [jj Q). Much 
attention has been attracted to studies of biomembranes 
with ion pumps and molecular motors, where membrane 
instabilities and synchronization effects have been ana- 
lyzed @, H, [I]- Here, a different class of distributed ac- 
tive molecular systems — formed by enzymes — is con- 
sidered. The catalytic activity of an allosteric enzyme 
protein is activated or inhibited by binding of small reg- 
ulatory molecules; the role of such regulatory molecules 
can be played by products of the same reaction [ld| • Pre- 
vious investigations of simple product-regulated enzymic 
systems [111 . Il2fl and enzymic networks [131 ] in small spa- 
tial volume with full diffusional mixing have shown that 
spontaneous synchronization of molecular turnover cycles 
can take place there. External molecular synchroniza- 
tion of enzymes of the photosensitive P-450 dependent 
monooxygenase system by periodic optical forcing has 
been experimentally demonstrated [l4j . 

In this Letter, spatiotemporal pattern formation in en- 
zymic arrays is investigated. In such systems, immobile 
enzymes are attached to a solid planar support immersed 
into a solution through which fresh substrate is supplied 
and product molecules arc continuously removed. Prod- 
uct molecules released by an enzyme diffuse through the 
solution and activate catalytic turnover cycles of neigh- 
bouring enzymes in the array. 

A simple stochastic model [H| of an enzyme as a cyclic 
machine (a stochastic phase oscillator), shown in Fig. 1, 
is used. Binding of a substrate molecule to an enzyme i 
initiates an ordered internal conformational motion, de- 
scribed by the conformational phase coordinate fa. The 
initial state corresponds to the phase fa = 0. The cat- 
alytic conversion event takes place and the product is 
released at the state fa inside the cycle. After that, the 
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FIG. 1: (Color online) A sketch of the model. 



conformational motion continues until the equilibrium 
state of the enzyme (fa = 1) is finally reached. Initi- 
ation of a turnover cycle is a random event, occurring 
at a certain probability rate. We assume that substrate 
is present in abundance, and its concentration is not af- 
fected by the reactions. Conformational motion inside 
the cycle is modeled as a stochastic diffusional drift pro- 
cess, described by equation fa = v + rji(t), where v is the 
mean drift velocity and rji(t) is an internal white noise 
with {rii(t)r)j(t')) — 2a5ijS(t — t') where a specifies inten- 
sity of intramolecular fluctuations. 

Allosterically activated enzymes possess a site on their 
surface where regulatory molecules can become bound. 
Binding of a regulatory molecule leads to conformational 
change that enhances catalytic activity of the enzyme. A 
regulatory molecule binds to an enzyme with rate con- 
stant [3 and dissociate from it with rate constant k. Bind- 
ing of a regulatory molecule at an enzyme raises its prob- 
ability to start a cycle from ao to a±. We assume that 
a regulatory molecule can bind to an enzyme only in its 
rest state and this molecule is released when the cycle 
is started. The role of regulatory molecules is played 
by product molecules of the same reaction. Immobile 
enzymes are randomly distributed in space with concen- 
tration c. Product diffuses at diffusion constant D and 
undergoes decay at rate constant 7. The characteristic 
diffusion length of product molecules is Idijf — yD/l- 

In our stochastic 2D simulations, the medium was dis- 
cretized into spatial cells (up to 256 x 256), each con- 
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FIG. 2: Stochastic (a,b) and mean-field (c,d) simulations of 
2D wave patterns; (a) t p = 0.14, c = 1, and j3 = 300, (b) 
t p = 0.25, c = 10, and /3 = 10, (c) t p = 0.14, c = 1, and 
(3 = 300, (d) t p = 0.34, c = 100, and /3 = 1.42. Other 
parameters are ao = 1, ai = 1000, « = 10, 7 = 10, cr = 0, 
D = 100. The linear size of the shown area is L = 40 Idiff in 
all panels. 



taining a number of enzyme molecules. The cells were 
so small that diffusional mixing of product molecules in 
a cell within the shortest characteristic time of the reac- 
tion could always take place. Each enzyme was described 
by the stochastic model given above; diffusion of product 
molecules was modeled as a random walk over a discrete 
cell lattice. The mean cycle time r = 1/v was chosen as 
the time unit (r = 1). Systems including up to 655 360 
enzymes were used in the simulations. 

Figure 2a, b (see also Videos 1 and 2 in ref. [15j) shows 
two typical examples of stochastic 2D simulations. Here, 
spatial distributions of product molecules are displayed. 
Waves of product concentration are propagating through 
the medium. In a peak of a wave, many locally present 
enzymes are simultaneously releasing product molecules. 
Since product release can take place only at a certain 
stage inside the cycle, this means that the cycles of en- 
zymes are locally synchronized. Not only regular wave 
structures, such as rotating spiral waves or target pat- 
terns (Fig. 2a), but also complex regimes of wave turbu- 
lence (Fig. 2b) have been observed. 

To understand and interpret stochastic simulation re- 
sults, an analytical study of the system in the mean-field 
approximation, which holds in the limit of high enzyme 
concentrations, has been performed. In this approxima- 
tion, the system is characterised by three continuous vari- 
ables no(r,t), ni(r,t) and m(r, t) which represent local 
concentrations of enzymes in the rest state without or 
with regulatory molecules attached (no and m) and local 
concentration of the product (m). For simplicity, internal 
fluctuations in enzymes are neglected (a = 0). Thus, all 
enzymes which have started their cycles at some time t 
would release their products at a definite time t+r p (with 
t p = 4>pl v ) an d finish their cycles, returning to the rest 
state, at time t + r. Therefore, the system is described 
by a set of three reaction-diffusion equations with time 



delays, 

-7j£- = pmno — Kfti — otirii (la) 

dno . . 

= -fjmriQ + nni - a no + a a n {t - t) 

+a 1 ni(t — t) (lb) 

dm 

= — fjmriQ + mix + ami - jm + aon (i — t p ) 
+aini(t — t p ) + DV 2 m. (lc) 

The system always has a uniform stationary state 
with certain concentrations no, n~i and m, which can be 
found as solutions of the respective algebraic equations. 
This state corresponds to the absence of synchronization. 
However, it may become unstable if allosteric activation 
is strong enough. To analyze stability, small perturba- 
tions dno, Sni and dm are added to the stationary state, 
equations (1) are linearized and their solutions are sought 
as Sno ~ Sni ~ 5m ~ exp (X q t — iqx) with X q = [i q + iuj q . 
Thus, each spatial mode with wavevector q is character- 
ized by its frequency u q and its rate of growth u q . The 
properties fx q and oj q are given by the roots of a charac- 
teristic equation which is determined by the linearization 
matrix of equations (1). The steady state becomes unsta- 
ble when at least one spatial mode with some wavenum- 
ber go starts to grow (fj, qo > 0). 

As the bifurcation parameter, coefficient f3 can be cho- 
sen. If regulatory molecules cannot bind to enzymes 
((3 = 0), feedback is absent and instabilities are not pos- 
sible. On the other hand, allosteric activation becomes 
strong if regulatory molecules can easily bind and, in this 
case, emergence of oscillations and wave patterns can be 
expected. Our bifurcation analysis reveals that, depend- 
ing on the parameters of the system, it can exhibit ei- 
ther a Hopf or a wave bifurcation [l6(. As a result of 
the Hopf bifurcation, uniform oscillations with q = de- 
velop. Because of the presence of delays in equations (1), 
the characteristic equation is nonpolynomial in terms of 
A and, generally, a number of oscillatory solutions with 
different frequencies u> are possible. Physically, such so- 
lutions correspond to formation of several synchronous 
enzymic groups. This effect has been previously exten- 
sively investigated for similar systems in small spatial 
volumes with full diffusional mixing [1 ll | and we shall 
not further discuss it here. The most robust uniform 
oscillations, which we consider, are characterized by the 
frequency w « 2it/t and correspond to the single-group 
synchronization. As the result of a wave bifurcation (also 
known as the Hopf bifurcation with a finite wave number 
[13]), the first unstable modes are traveling waves with 
a certain wavenumber go- Figure 3 shows the bifurca- 
tion diagram in the parameter plane (r p , j3). Note the 
presence of a codimension-2 bifurcation point where the 
boundaries of the Hopf and the wave bifurcations join. 

To investigate nonlinear dynamics of the system, nu- 
merical simulations of equations (1) have been performed 
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FIG. 3: Phase diagram (q = 1, ai = 1000, k = 10, 7 = 10, 
c = 100, D = 1000). The Hopf bifurcation (solid line) and the 
wave bifurcation (dash-dotted line) boundaries are displayed. 
Gray lines show instability of the stationary state with respect 
to development of uniform oscillations with two (dashed) and 
three (dotted) groups in the well-mixed case. Lines separating 
parameter domains with different kinds of patterns are hand- 
drawn, based on numerical simulations. 



FIG. 4: Spatiotemporal patterns in a ID system (in each 
panel, the vertical axis is time, running down, and the hor- 
izontal axis is the coordinate). The upper two rows are 
stochastic simulations (a — 0) with concentrations c = 1 and 
c = 10, the bottom row shows mean-field simulations with 
c = 100. (a) t p = 0.3, (3 = 95/c, (b) t p = 0.14, (3 = 260/c, (c) 
t p = 0.22, (3 = 600/c, and (d) t p = 0.16, (3 = 300/c. Other 



parameters as in Fig. 
Idiff- 



3; the system size shown is L = 51 



[16j . The explicit Euler integration method has been 
used; no-flux boundary conditions were applied. Results 
of ID simulations are summarized in Fig. 3 and examples 
of typical observed patterns are shown in Fig. 4. Stand- 
ing waves (Fig. 4a) develop when the boundary of the 
wave bifurcation (dash-dotted curve) is crossed and uni- 
form oscillations are observed above the boundary of the 
Hopf bifurcation. Near the codimension-2 point, more 
complex behavior was found. This included rippled os- 
cillations (Fig. 4b), self-organized pacemakers (Fig. 4c) 
and modulated traveling waves (Fig. 4d). The observed 
patterns are similar to those previously found in reaction- 
diffusion systems with the wave bifurcation (l8| . In the 
right upper corner of the diagram in Fig. 3, higher fre- 
quency oscillations with several synchronous groups take 
place. 

Two-dimensional simulations of reaction-diffusion 
equations (1) with time delay have been performed for 
selected parameter values. In 2D simulations, sponta- 
neously developing concentric waves (target patterns) 
and spiral waves have been observed; target patterns 
were however unstable and evolved into pairs of rotat- 
ing spiral waves (Fig. 2c and Video 3 [15(). Complex 
wave regimes, which can be qualitatively characterized 
as turbulence of standing waves, have also been observed 
(Fig. 2d and Video 4 [Tj|). 

The mean-field approximation is based on neglect- 
ing statistical fluctuations in concentrations of reacting 
species ll| and, therefore, it should hold in the high 
concentration limit. In Fig. 4, two upper panel rows 
display spatiotemporal patterns which are observed in 



stochastic simulations with parameter values correspond- 
ing to the respective mean-field simulations. To compare 
mean-field simulations with different enzyme densities, 
the following property of equations (1) can be used: in- 
troducing relative concentrations uq — tiq/c, ri\ = n\jc 
and to = m/c, it can be noticed that they obey the same 
equations, but with a rescaled coefficient (3 = (3c. Thus, 
essentially the same patterns are observed as long as 
the parameter combination (3c remains constant. In the 
stochastic simulations in Fig. 4, the coefficient (3 has been 
increased to compensate for a decrease in the enzyme 
concentration. For larger enzyme concentrations, good 
agreement between mean-field predictions and stochas- 
tic simulations has been found. In the mean-field equa- 
tions (1), intramolecular fluctuations are not taken into 
account {a = and therefore each turnover cycle has 
the same fixed duration r). Stochastic simulations have 
been, however, also performed when such fluctuations 
were present. Synchronization waves could still be found 
even at internal noise levels which corresponded to the 
mean relative dispersion £ of turnover times of about 10% 

(with £ = (St 2 ) 1 ' 2 It ~ (2ctt) 1/2 ). 

Although the emphasis in this Letter is on the phenom- 
ena in two-dimensional enzymic arrays, analogous effects 
should be expected for three-dimensional systems repre- 
senting aqueous enzymic solutions. The linear stability 
analysis, yielding Hopf and wave bifurcation boundaries 
(see Fig. 3), is valid also for the 3D geometry. We have 
performed preliminary stochastic simulations for thin so- 
lution layers with high enzyme concentrations and could 
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observe synchronization patterns similar to those found 
for the enzymic arrays. 

A product molecule, released by an enzyme, diffuses 
in the solution until it either binds, as a regulatory 
molecule, to another enzyme or undergoes a decay. Here, 
it should be taken into account that a regulatory molecule 
can bind to an allosteric enzyme only at a certain bind- 
ing site of characteristic radius R. Using the theory of 
diffusion-controlled reactions, the average time t tra nsit 
after which a regulatory product would find a binding 
site of one of the enzymes can be roughly estimated 
ll| as ttransit = X/cDR, if enzymes are uniformly dis- 
tributed inside the reaction volume with concentration 
c. Therefore, binding typically occurs within the dis- 
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from the point 



where a molecule is released. Obviously, it can only 
take place if the product molecule has not undergone 
decay until that moment, i.e. if ^ttransit < 1- This 
condition puts a restriction on the enzyme concentration 
c, which must be higher than the critical concentration 
c* = j/DR. Choosing 7 = 10 3 s -1 , D = 1(T 5 cnr^s" 1 
and R = 10~ 7 cm, the critical enzyme concentration is 
c* = 10 15 cm~ 3 = fCP 6 M. A similar estimate can be 
obtained when enzymes are immobilized on a plane im- 
mersed into a reactive solution; in this case the mean dis- 
tance between the enzymes on the plane should be less 



than l c = {Rldiff 



,1/2 



Although the required en- 



zyme concentrations are relatively large, they are within 
the range characteristic for biological cells (glycolytic en- 
zymes are present [l9j in a cell at even higher concentra- 
tion of more than 10 -5 M). The characteristic temporal 
period of developing patterns is determined by the en- 
zyme turnover time r, which typically varies from mil- 
liseconds to seconds. The characteristic length scale of 
developing wave patterns is determined by the diffusion 
length Idiff, which can vary under these conditions from 
a fraction of a micrometer to tens of micrometers. 

Our analysis shows that spontaneous molecular syn- 
chronization of allosteric product-activated enzymes can 
be observed in enzymic arrays. Artificial arrays formed 
by immobilized protein machines (molecular motors) are 
already used in experiments on active nanoscale trans- 
port (see [20]). Many enzymes in biological cells are 
membrane-bound, thus forming natural enzymic arrays. 
Similar phenomena are possible in dense enzyme solu- 



tions. In the study by Petty et al. [2l[, traveling waves 
of NAD(P)H and proton concentrations with the wave- 
length of about a micrometer were observed inside neu- 
trophil cells. These metabolic waves had the temporal 
period of about 300 ms, which is by two orders of magni- 
tude shorter than the characteristic period of glycolytic 
oscillations in the cells and lies closer to the time scales 
of turnover cycles of individual enzymes. An intriguing 
question, requiring further detailed analysis, is whether 
molecular synchronization waves may have already been 



seen in these experiments. 

Molecular synchronization waves are principally dif- 
ferent from classical concentration waves in reaction- 
diffusion systems. Under synchronization conditions, 
internal conformational states of individual enzyme 
molecules in their turnover cycles become strongly cor- 
related. In optics, a similar situation is found when a 
transition to coherent laser generation has taken place. 
Our theoretical analysis may open a way to the investiga- 
tions of a new class of spatio-temporal pattern formation 
in chemically active molecular systems. 
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